DEVICES = [
    'beaglebone-fan',
    'beaglebone-compressor',
    'beaglebone-pump',
    'beaglebone-refrigerator',
    'mafaulda-a',
    'mafaulda-b'
]

DEVICE = DEVICES[5]
T_WAVEFORM = 1  # (1 = MaufaulDa, 5 = others)
T_SEC = T_WAVEFORM
NFFT = 2**14  # (2**10 = MaufaulDa, 2**14 = others)
F_LIMIT = 3000 # (3000 = MaufaulDa, None = others)

import os
import pandas as pd
import numpy as np
import matplotlib.pylab as plt

from tabulate import tabulate
from IPython.display import Markdown, HTML
from tqdm.notebook import tqdm

from typing import List, Tuple
from scipy.signal import find_peaks, butter, lfilter
from tsfel.feature_extraction.features import fundamental_frequency

from zipfile import ZipFile
import sys
sys.path.append('../')
from vibrodiagnostics import mafaulda, selection, discovery, models


def beaglebone_measurement(filename: str, fs: int) -> Tuple[str, pd.DataFrame]:
    g = 9.81
    milivolts = 1800
    resolution = 2**12
    columns = ['x', 'y', 'z']
    ts = pd.read_csv(filename, delimiter='\t', index_col=False, header=None, names=columns)
        
    # Calculate amplitude in m/s^2 Beaglebone Black ADC and ADXL335 resolution (VIN 1.8V, 12bits)
    for dim in columns:
        ts[dim] = ts[dim] * (milivolts / resolution)  # ADC to mV
        ts[dim] = (ts[dim] / 180) * g                 # mV to m/s^2 (180 mV/g)
        ts[dim] -= ts[dim].mean()

    ts['t'] = ts.index * (1 / fs)
    ts.set_index('t', inplace=True)
    return (os.path.basename(filename), ts, fs, ts.columns)  # last is feature columns


def beaglebone_dataset(filenames: List[str], fs: int) -> List[Tuple[str, pd.DataFrame]]:
    dataset = []
    for filename in filenames:
        name, ts, fs, cols = beaglebone_measurement(filename, fs)
        dataset.append((name, ts))
    return dataset


def lowpass_filter(data, cutpoint, fs, order=5):
    b, a = butter(order, cutpoint, fs=fs, btype='lowpass')
    y = lfilter(b, a, data)
    return y


def mafaulda_dataset(
        place,
        features_path =  '../../datasets/features_data/',
        mafaulda_path='../../datasets/MAFAULDA.zip',
        rpm=2500,
        lowpass_hz=10000):

    metadata_filename = os.path.join(features_path, selection.MAFAULDA_METADATA)
    faults = {
        'A': {
            'normal': 'normal',
            'imbalance': 'imbalance',
            'horizontal-misalignment': 'misalignment',
            'vertical-misalignment': 'misalignment',
            'underhang-outer_race': 'outer race fault',
            'underhang-cage_fault': 'cage fault',
            'underhang-ball_fault': 'ball fault'
        },
        'B': {
            'normal': 'normal',
            'imbalance': 'imbalance',
            'horizontal-misalignment': 'misalignment',
            'vertical-misalignment': 'misalignment',
            'overhang-cage_fault': 'cage fault',
            'overhang-ball_fault': 'ball fault',
            'overhang-outer_race': 'outer race fault'
        }
    }
    bearings = {
        'A': ['ax', 'ay', 'az'],
        'B': ['bx', 'by', 'bz']
    }

    metadata = pd.read_csv(metadata_filename, index_col='filename')
    metadata.reset_index(inplace=True)
    metadata = metadata[metadata['fault'].isin(tuple(faults[place]))]
    metadata = models.fault_labeling(metadata, faults[place])
    files = pd.DataFrame()
    # Worst severity and mid rpm
    for name, group in metadata[(metadata['rpm'] >= rpm)].groupby(by='fault', observed=False):
        files = pd.concat([
            files,
            group[
                group['severity_level'] == group['severity_level'].max()
            ].sort_values(by='rpm', ascending=True).head(1)
        ])
    ordering = {
        'normal': 0,
        'misalignment': 1,
        'imbalance': 2,
        'cage fault': 3,
        'ball fault': 4,
        'outer race fault': 5,
    }
    source = ZipFile(mafaulda_path)
    dataset = len(files) * [0]
    for index, file in files.iterrows():
        ts = mafaulda.csv_import(source, file['filename'])
        ts = ts[bearings[place]]
        ts.columns = ts.columns.str.extract(r'(\w)$')[0]
        for axis in ts.columns:
            ts[axis] = lowpass_filter(ts[axis], lowpass_hz, file['fs'])
        pos = ordering[file['fault']]
        dataset[pos] = ((file['fault'] + ' (' + file['filename'] +')', ts))

    return dataset

if DEVICE == 'beaglebone-fan':
    Fs = 2500
    path = '../../inspections/fan/'
    files = [
        '1_still.tsv', '2_still.tsv', '3_still.tsv',
        '1_up.tsv', '2_up.tsv', '3_up.tsv',
        '1_down.tsv', '2_down.tsv', '3_down.tsv'
    ]
    files = [os.path.join(path, name) for name in files]
    DATASET = beaglebone_dataset(files, Fs)

elif DEVICE == 'beaglebone-compressor':
    Fs = 2500
    path = '../../inspections/datacentres/shc3/'
    files = [
        'k3_1.tsv', 'k3_2.tsv', 'k3_3.tsv', 'k3_4.tsv',
        'k5_1.tsv', 'k5_2.tsv', 'k5_3.tsv', 'k5_4.tsv'
    ]
    files = [os.path.join(path, name) for name in files]
    DATASET = beaglebone_dataset(files, Fs)

elif DEVICE == 'beaglebone-pump':
    Fs = 2500
    path = '../../inspections/bvs/'
    files = [
        'bvs_1_hore.tsv', 'bvs_2_hore.tsv' 
        #, 'bvs_3_motor.tsv', 'bvs_4_motor.tsv'
    ]    
    files = [os.path.join(path, name) for name in files]
    DATASET = beaglebone_dataset(files, Fs)

elif DEVICE == 'beaglebone-refrigerator':
    Fs = 2500
    path = '../../inspections/home-refrigerator/'
    files = [
        'ch1.tsv', 'ch2.tsv', 'ch3.tsv', 'ch4.tsv', 'ch5.tsv'
    ]    
    files = [os.path.join(path, name) for name in files]
    DATASET = beaglebone_dataset(files, Fs)

elif DEVICE == 'mafaulda-a':
    Fs = mafaulda.FS_HZ
    DATASET = mafaulda_dataset(place='A')

elif DEVICE == 'mafaulda-b':
    Fs = mafaulda.FS_HZ
    DATASET = mafaulda_dataset(place='B')

for name, ts in DATASET:
    display(Markdown(f'**{name}**'))
    ts.info()
    print()

<class 'pandas.core.frame.DataFrame'>
Index: 250000 entries, 0.0 to 4.999980000000001
Data columns (total 3 columns):
 #   Column  Non-Null Count   Dtype  
---  ------  --------------   -----  
 0   x       250000 non-null  float64
 1   y       250000 non-null  float64
 2   z       250000 non-null  float64
dtypes: float64(3)
memory usage: 7.6 MB

<class 'pandas.core.frame.DataFrame'>
Index: 250000 entries, 0.0 to 4.999980000000001
Data columns (total 3 columns):
 #   Column  Non-Null Count   Dtype  
---  ------  --------------   -----  
 0   x       250000 non-null  float64
 1   y       250000 non-null  float64
 2   z       250000 non-null  float64
dtypes: float64(3)
memory usage: 7.6 MB

<class 'pandas.core.frame.DataFrame'>
Index: 250000 entries, 0.0 to 4.999980000000001
Data columns (total 3 columns):
 #   Column  Non-Null Count   Dtype  
---  ------  --------------   -----  
 0   x       250000 non-null  float64
 1   y       250000 non-null  float64
 2   z       250000 non-null  float64
dtypes: float64(3)
memory usage: 7.6 MB

<class 'pandas.core.frame.DataFrame'>
Index: 250000 entries, 0.0 to 4.999980000000001
Data columns (total 3 columns):
 #   Column  Non-Null Count   Dtype  
---  ------  --------------   -----  
 0   x       250000 non-null  float64
 1   y       250000 non-null  float64
 2   z       250000 non-null  float64
dtypes: float64(3)
memory usage: 7.6 MB

<class 'pandas.core.frame.DataFrame'>
Index: 250000 entries, 0.0 to 4.999980000000001
Data columns (total 3 columns):
 #   Column  Non-Null Count   Dtype  
---  ------  --------------   -----  
 0   x       250000 non-null  float64
 1   y       250000 non-null  float64
 2   z       250000 non-null  float64
dtypes: float64(3)
memory usage: 7.6 MB

<class 'pandas.core.frame.DataFrame'>
Index: 250000 entries, 0.0 to 4.999980000000001
Data columns (total 3 columns):
 #   Column  Non-Null Count   Dtype  
---  ------  --------------   -----  
 0   x       250000 non-null  float64
 1   y       250000 non-null  float64
 2   z       250000 non-null  float64
dtypes: float64(3)
memory usage: 7.6 MB

for name, ts in DATASET:
    display(Markdown(f'**{name}**'))
    display(tabulate(ts.describe(), headers='keys', tablefmt='html'))
    ts.boxplot(grid=True)
    plt.show()

from statsmodels.tsa.stattools import adfuller
from statsmodels.api import qqplot
from scipy.stats import kstest

normality_tests = []
for name, ts in DATASET:
    for x in ts.columns:
        p_value = kstest(ts[x], 'norm').pvalue
        test = {'name': name, 'axis': x, 'p-value': p_value, 'not-normal': p_value < 0.05}
        normality_tests.append(test)

normality_tests = pd.DataFrame.from_records(normality_tests)
print(normality_tests.value_counts('not-normal'))
normality_tests.describe()

not-normal
True    18
Name: count, dtype: int64

name, ts = DATASET[0]
fig, ax = plt.subplots(1, len(ts.columns), figsize=(10, 4))
for i, x in enumerate(ts.columns):
    qqplot(ts[x], line='45', ax=ax[i], marker='.', alpha=0.5)
    ax[i].set_title(f'Axis: {x}')

plt.tight_layout()
print(name)
plt.show()

normal (normal/43.6224.csv)

stationarity_tests = []
for name, ts in tqdm(DATASET):
    for x in ts.columns:
        result = adfuller(ts[x].loc[T_WAVEFORM:T_WAVEFORM+1])
        p_value = result[1]
        test = {
            'name': name,
            'axis': x,
            'statistic': result[0],
            'p-value': p_value,
            'stationary': p_value < 0.001
        }
        stationarity_tests.append(test)

stationarity_tests = pd.DataFrame.from_records(stationarity_tests)
print(stationarity_tests.value_counts('stationary'))
stationarity_tests['p-value'].describe()

  0%|          | 0/6 [00:00<?, ?it/s]

stationary
True    18
Name: count, dtype: int64

count    1.800000e+01
mean     5.255316e-15
std      1.883517e-14
min      0.000000e+00
25%      1.288097e-29
50%      2.940757e-23
75%      2.062478e-18
max      7.956672e-14
Name: p-value, dtype: float64

name, ts = DATASET[0]
fig, ax = plt.subplots(1, len(ts.columns), figsize=(10, 4))
for i, x in enumerate(ts.columns):
    ax[i].acorr(ts[x], maxlags=50)
    ax[i].set_title(f'Axis: {x}')

plt.tight_layout()
print(name)
plt.show()

normal (normal/43.6224.csv)

for name, ts in DATASET:
    display(Markdown(f'**{name}**'))
    axis = ts.columns
    ax = ts[axis].hist(figsize=(15, 3), grid=True, bins=100, layout=(1, 3), edgecolor='black', linewidth=0.5)
    plt.show()

for name, ts in DATASET:
    display(Markdown(f'**{name}**'))
    axis = ts.columns
    
    ax = ts[axis].plot(figsize=(20, 8), grid=True, subplots=True)
    for i, axname in enumerate(axis):
        ax[i].set_xlabel('Time [s]')
        ax[i].set_ylabel(f'Amplitude ({axname}) [m/s^2]')
    plt.show()               # plt.savefig('waveform.png')

for name, ts in DATASET:
    axis = ts.columns
    display(Markdown(f'**{name}**'))
    ax = (ts[axis].iloc[int(T_WAVEFORM*Fs):int(T_WAVEFORM*Fs)+Fs]
                  .plot(figsize=(20, 10), grid=True, subplots=True))
    
    for i, axname in enumerate(axis):
        ax[i].set_xlabel('Time [s]')
        ax[i].set_ylabel(f'Amplitude ({axname}) [m/s^2]')
        plt.show()      # plt.savefig('waveform_zoom.png')

fig, ax = plt.subplots(len(DATASET), 3, figsize=(12, 15), sharex=True)

for idx, df in enumerate(DATASET):
    name, ts = df
    columns = ts.columns
    ax[idx][1].set_title(name)
    ax[idx][0].set_ylabel('Amplitude [m/s^2]')

    for pos, axis in enumerate(columns):
        data = ts[axis].loc[T_WAVEFORM:T_WAVEFORM+0.3]
        ax[idx][pos].plot(data.index, data, linewidth=1, color='darkblue')
        ax[idx][pos].grid()
    

plt.tight_layout()
plt.show()

def spectogram(x, debug=True):
    fig, ax = plt.subplots(figsize=(15, 4))
    cmap = plt.get_cmap('inferno')
    pxx, freqs, t, im = plt.specgram(
        x, NFFT=NFFT, Fs=Fs,
        detrend='mean',
        mode='magnitude', scale='dB',
        cmap=cmap, vmin=-60
    )
    fig.colorbar(im, aspect=20, pad=0.04)
    ax.set_xlabel('Time [s]')
    ax.set_ylabel('Frequency [Hz]')
    mafaulda.resolution_calc(Fs, NFFT)
    return freqs, pxx


def window_idx(t):
    return (Fs * t) // NFFT + 1


def spectrum_slice(freqs, Pxx, t):
    fig, ax = plt.subplots(2, 1, figsize=(20, 8))
    n = window_idx(t)

    dB = 20 * np.log10(Pxx.T[n] / 0.000001)
    ax[0].plot(freqs, dB)      # 1 dB = 1 um/s^2
    ax[0].grid(True)
    ax[0].set_xlabel('Frequency [Hz]')
    ax[0].set_ylabel('Amplitude [dB]')
    
    ax[1].plot(freqs, Pxx.T[n])
    ax[1].grid(True)
    ax[1].set_xlabel('Frequency [Hz]')
    ax[1].set_ylabel('Amplitude [m/s^2]')
    return n


def get_max_frequency(freqs, Pxx, i):
    max_freq = freqs[np.argmax(Pxx.T[i])]
    return max_freq


def get_peaks(freqs, Pxx, i, top=5):
    amplitudes = Pxx.T[i]
    peaks, _ = find_peaks(amplitudes, distance=3)

    fundamental = get_max_frequency(freqs, Pxx, i)
    f_top = freqs[peaks[np.argsort(amplitudes[peaks])]][::-top]
    y_top = np.sort(amplitudes[peaks])[::-top]

    return pd.DataFrame({
        'f': f_top,
        'y': y_top,
        '1x': f_top / fundamental 
    })


def butter_bandpass_filter(data, lowcut, highcut, fs, order=5):
    b, a = butter(order, [lowcut, highcut], fs=fs, btype='band')
    y = lfilter(b, a, data)
    return y


def get_spectrograms(DATASET: List[pd.DataFrame], axis: str) -> list:
    spectrograms = []

    for name, ts in DATASET:
        base_freq = fundamental_frequency(ts[axis], Fs)
        display(Markdown(f'**{name}** *({axis.upper()} axis, Fundamental = {base_freq:.4f} Hz)*'))
        
        freqs, Pxx = spectogram(ts[axis])
        spectrograms.append((name, freqs, Pxx))
        plt.show()          # plt.savefig(f'x_axis_fft_{NFFT}.png')
    
    return spectrograms


def show_spectrogram_detail(spectrograms: list, axis: str, t: float):
    for name, freqs, Pxx in spectrograms:
        display(Markdown(f'**{name}** ({axis.upper()} axis @ {t}s)'))
        i_window = spectrum_slice(freqs, Pxx, t)
        plt.show()           #plt.savefig(f'x_axis_fft_{NFFT}_at_{T_SEC}s.png')


def show_mms_peaks(spectrograms: list, axis: str, t: float):
    for name, freqs, Pxx in spectrograms:
        display(Markdown(f'**{name}** ({axis.upper()} axis @ {t}s)'))
    
        i_window = window_idx(t)
        peaks = discovery.mms_peak_finder(Pxx.T[i_window])
        
        fig, ax = plt.subplots(1, 1, figsize=(15, 3))
        ax.grid(True)
        ax.plot(freqs, Pxx.T[i_window])
        ax.scatter(freqs[peaks], Pxx.T[i_window][peaks], marker='^', color='red')
        ax.set_xlabel('Frequency [Hz]')
        
        plt.show()


def show_harmonic_series(spectrograms: list, axis: str, t: float):
    # https://stackoverflow.com/questions/1982770/changing-the-color-of-an-axis
    for name, freqs, Pxx in spectrograms:
        display(Markdown(f'**{name}** ({axis.upper()} axis @ {t}s)'))
    
        i_window = window_idx(t)
        h_series = discovery.harmonic_series_detection(freqs, Pxx.T[i_window], Fs, NFFT)
    
        # Find best (sum of harmonics' amplitudes in the largest)
        max_harmonic_amp_idx = np.argmax([
            sum([h[1] for h in s]) / len(s)
            for s in h_series
        ])
        best_harmonic_series = pd.DataFrame(
            h_series[max_harmonic_amp_idx],
            columns=['Frequency [Hz]', 'Amplitude [m/s^2]']
        )
        best_harmonic_series.index += 1
        display(tabulate(best_harmonic_series, headers='keys', tablefmt='html'))
    
        # Plot found harmonic series
        fig, ax = plt.subplots(1, 8, figsize=(30, 4))
        for i in range(8):
            s = h_series[i+1]
            if i == max_harmonic_amp_idx:
                ax[i].xaxis.label.set_color('red')
    
            ax[i].plot(freqs, Pxx.T[i_window])
            ax[i].scatter([x[0] for x in s], [x[1] for x in s], marker='^', color='red')
            ax[i].set_xlabel('Frequency [Hz]')
    
        plt.show()

def show_spectra_largest_amplitudes(spectrograms: list, axis: str, t: float):
    for name, freqs, Pxx in spectrograms:
        display(Markdown(f'**{name}** ({axis.upper()} axis @ {t}s)'))

        i_window = window_idx(t)
        x_fundamental = get_max_frequency(freqs, Pxx, i_window)
        peaks = get_peaks(freqs, Pxx, i_window)
        
        display(Markdown(f'- *Fundamental frequency:* {x_fundamental} Hz'))
        display(tabulate(peaks.head(5), headers='keys', tablefmt='html'))


def compare_limited_specrograms(spectrograms: list, axis: str, t: float):
    fig, ax = plt.subplots(len(DATASET), 1, figsize=(20, 20), sharey=True)
    i = 0
    for name, ts in DATASET:
        signal = ts[axis].loc[t:t+NFFT/Fs].to_numpy()
        pxx = np.abs(np.fft.rfft(signal) / len(signal)) 
        freqs = np.fft.fftfreq(len(signal), d=1/Fs)[:len(pxx)]
        #ilast = len(freqs[freqs < F_LIMIT])
        
        ax[i].plot(freqs, pxx)
        ax[i].grid(True)
        ax[i].set_xlabel('Frequency [Hz]')
        ax[i].set_ylabel('Amplitude [m/s^2]')
        ax[i].set_xlim(0, F_LIMIT)
        ax[i].set_ylim(0, 0.4)
        ax[i].set_title(name)
        i += 1


def spectrogram_energy_left_cumulative(spectrograms: list, axis: str, t: float):
    fig, ax = plt.subplots(len(DATASET), 1, figsize=(20, 20), sharey=True)
    i = 0
    for name, ts in DATASET:
        signal = ts[axis].loc[t:t+NFFT/Fs].to_numpy()
        pxx = np.abs(np.fft.rfft(signal) / len(signal)) 
        freqs = np.fft.fftfreq(len(signal), d=1/Fs)[:len(pxx)]
        
        ax[i].plot(freqs, np.cumsum(pxx) / np.sum(pxx))
        ax[i].grid(True)
        ax[i].set_xlabel('Frequency [Hz]')
        ax[i].set_ylabel('Cumulative energy [%]')
        #ax[i].set_xlim(0, 10000)
        ax[i].set_title(name)
        i += 1

compare_limited_specrograms(DATASET, 'x', T_SEC)
plt.tight_layout()
plt.show()

compare_limited_specrograms(DATASET, 'y', T_SEC)
plt.tight_layout()
plt.show()

compare_limited_specrograms(DATASET, 'z', T_SEC)
plt.tight_layout()
plt.show()

spectrogram_energy_left_cumulative(DATASET, 'x', T_SEC)
plt.tight_layout()
plt.show()

spectrogram_energy_left_cumulative(DATASET, 'y', T_SEC)
plt.tight_layout()
plt.show()

spectrogram_energy_left_cumulative(DATASET, 'z', T_SEC)
plt.tight_layout()
plt.show()

x_spectra = get_spectrograms(DATASET, 'x')

Window size: 16384
Heinsenberg box
	Time step: 327.68 ms
	Frequency step: 3.0517578125 Hz

Window size: 16384
Heinsenberg box
	Time step: 327.68 ms
	Frequency step: 3.0517578125 Hz

Window size: 16384
Heinsenberg box
	Time step: 327.68 ms
	Frequency step: 3.0517578125 Hz

Window size: 16384
Heinsenberg box
	Time step: 327.68 ms
	Frequency step: 3.0517578125 Hz

Window size: 16384
Heinsenberg box
	Time step: 327.68 ms
	Frequency step: 3.0517578125 Hz

Window size: 16384
Heinsenberg box
	Time step: 327.68 ms
	Frequency step: 3.0517578125 Hz

show_spectrogram_detail(x_spectra, 'x', T_SEC)

show_mms_peaks(x_spectra, 'x', T_SEC)

# show_harmonic_series(x_spectra, 'x', T_SEC)

show_spectra_largest_amplitudes(x_spectra, 'x', T_SEC)

y_spectra = get_spectrograms(DATASET, 'y')

Window size: 16384
Heinsenberg box
	Time step: 327.68 ms
	Frequency step: 3.0517578125 Hz

Window size: 16384
Heinsenberg box
	Time step: 327.68 ms
	Frequency step: 3.0517578125 Hz

Window size: 16384
Heinsenberg box
	Time step: 327.68 ms
	Frequency step: 3.0517578125 Hz

Window size: 16384
Heinsenberg box
	Time step: 327.68 ms
	Frequency step: 3.0517578125 Hz

Window size: 16384
Heinsenberg box
	Time step: 327.68 ms
	Frequency step: 3.0517578125 Hz

Window size: 16384
Heinsenberg box
	Time step: 327.68 ms
	Frequency step: 3.0517578125 Hz

show_spectrogram_detail(y_spectra, 'y', T_SEC)

show_mms_peaks(y_spectra, 'y', T_SEC)

# show_harmonic_series(y_spectra, 'y', T_SEC)

show_spectra_largest_amplitudes(y_spectra, 'y', T_SEC)

z_spectra = get_spectrograms(DATASET, 'z')

Window size: 16384
Heinsenberg box
	Time step: 327.68 ms
	Frequency step: 3.0517578125 Hz

Window size: 16384
Heinsenberg box
	Time step: 327.68 ms
	Frequency step: 3.0517578125 Hz

Window size: 16384
Heinsenberg box
	Time step: 327.68 ms
	Frequency step: 3.0517578125 Hz

Window size: 16384
Heinsenberg box
	Time step: 327.68 ms
	Frequency step: 3.0517578125 Hz

Window size: 16384
Heinsenberg box
	Time step: 327.68 ms
	Frequency step: 3.0517578125 Hz

Window size: 16384
Heinsenberg box
	Time step: 327.68 ms
	Frequency step: 3.0517578125 Hz

show_spectrogram_detail(z_spectra, 'z', T_SEC)

show_mms_peaks(z_spectra, 'z', T_SEC)

# show_harmonic_series(z_spectra, 'z', T_SEC)

show_spectra_largest_amplitudes(z_spectra, 'z', T_SEC)

axis = ['x', 'y', 'z']
for name, ts in DATASET:
    display(Markdown(f'**{name}**'))
    ts[axis].hist(figsize=(10, 5), grid=True, bins=50)
    plt.show()

axis = ['x', 'y', 'z']
for name, ts in DATASET:
    display(Markdown(f'**{name}**'))
    ts[axis].boxplot(figsize=(10, 5))
    plt.show()

for name, ts in DATASET:
    display(Markdown(f'**{name}**'))
    fig, ax = plt.subplots(1, 3, figsize=(20, 4))

    for i, col in enumerate([('x', 'y'), ('x', 'z'), ('y', 'z')]):
        ax[i].scatter(ts[col[0]], ts[col[1]], s=1)
        ax[i].grid(True)
        ax[i].set_xlabel(col[0].upper())
        ax[i].set_ylabel(col[1].upper())
        ax[i].grid(True)
    plt.show()       # plt.savefig('orbitals.png')

x_spectra_by_name = {spec[0]: spec for spec in x_spectra}
y_spectra_by_name = {spec[0]: spec for spec in y_spectra}
z_spectra_by_name = {spec[0]: spec for spec in z_spectra}
t = 5
space = 5

for name, ts in DATASET:
    display(Markdown(f'**{name}**'))
    fig, ax = plt.subplots(1, 3, figsize=(20, 4))

    name, freqs, Pxx = x_spectra_by_name[name]
    x_fundamental = get_max_frequency(freqs, Pxx, window_idx(t))
    name, freqs, Pxx = y_spectra_by_name[name]
    y_fundamental = get_max_frequency(freqs, Pxx, window_idx(t))

    name, freqs, Pxx = z_spectra_by_name[name]
    z_fundamental = get_max_frequency(freqs, Pxx, window_idx(t))

    try:
        ts['x_1x'] = butter_bandpass_filter(ts['x'], x_fundamental - space, x_fundamental + space, Fs)
        ts['y_1x'] = butter_bandpass_filter(ts['y'], y_fundamental - space, y_fundamental + space, Fs)
        ts['z_1x'] = butter_bandpass_filter(ts['z'], z_fundamental - space, z_fundamental + space, Fs)
    except ValueError:
        continue
    
    for i, col in enumerate([('x_1x', 'y_1x'), ('x_1x', 'z_1x'), ('y_1x', 'z_1x')]):
        ax[i].scatter(ts[col[0]], ts[col[1]], s=1)
        ax[i].grid(True)
        ax[i].set_xlabel(col[0].upper())
        ax[i].set_ylabel(col[1].upper())
        ax[i].grid(True)
    
    plt.show()       # plt.savefig('orbitals_1x.png')

---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
/home/miroslav/fiit-stu/Ing/3-semester/DP/masters-thesis/notebooks/ExploratoryAnalysis/EDA.ipynb Cell 63 line 1
      <a href='vscode-notebook-cell:/home/miroslav/fiit-stu/Ing/3-semester/DP/masters-thesis/notebooks/ExploratoryAnalysis/EDA.ipynb#Y114sZmlsZQ%3D%3D?line=8'>9</a> fig, ax = plt.subplots(1, 3, figsize=(20, 4))
     <a href='vscode-notebook-cell:/home/miroslav/fiit-stu/Ing/3-semester/DP/masters-thesis/notebooks/ExploratoryAnalysis/EDA.ipynb#Y114sZmlsZQ%3D%3D?line=10'>11</a> name, freqs, Pxx = x_spectra_by_name[name]
---> <a href='vscode-notebook-cell:/home/miroslav/fiit-stu/Ing/3-semester/DP/masters-thesis/notebooks/ExploratoryAnalysis/EDA.ipynb#Y114sZmlsZQ%3D%3D?line=11'>12</a> x_fundamental = get_max_frequency(freqs, Pxx, window_idx(t))
     <a href='vscode-notebook-cell:/home/miroslav/fiit-stu/Ing/3-semester/DP/masters-thesis/notebooks/ExploratoryAnalysis/EDA.ipynb#Y114sZmlsZQ%3D%3D?line=12'>13</a> name, freqs, Pxx = y_spectra_by_name[name]
     <a href='vscode-notebook-cell:/home/miroslav/fiit-stu/Ing/3-semester/DP/masters-thesis/notebooks/ExploratoryAnalysis/EDA.ipynb#Y114sZmlsZQ%3D%3D?line=13'>14</a> y_fundamental = get_max_frequency(freqs, Pxx, window_idx(t))

/home/miroslav/fiit-stu/Ing/3-semester/DP/masters-thesis/notebooks/ExploratoryAnalysis/EDA.ipynb Cell 63 line 3
     <a href='vscode-notebook-cell:/home/miroslav/fiit-stu/Ing/3-semester/DP/masters-thesis/notebooks/ExploratoryAnalysis/EDA.ipynb#Y114sZmlsZQ%3D%3D?line=37'>38</a> def get_max_frequency(freqs, Pxx, i):
---> <a href='vscode-notebook-cell:/home/miroslav/fiit-stu/Ing/3-semester/DP/masters-thesis/notebooks/ExploratoryAnalysis/EDA.ipynb#Y114sZmlsZQ%3D%3D?line=38'>39</a>     max_freq = freqs[np.argmax(Pxx.T[i])]
     <a href='vscode-notebook-cell:/home/miroslav/fiit-stu/Ing/3-semester/DP/masters-thesis/notebooks/ExploratoryAnalysis/EDA.ipynb#Y114sZmlsZQ%3D%3D?line=39'>40</a>     return max_freq

IndexError: index 16 is out of bounds for axis 0 with size 15

t = 5
space = 8

for name, ts in DATASET:
    display(Markdown(f'**{name}**'))

    name, freqs, Pxx = x_spectra_by_name[name]
    x_fundamental = get_max_frequency(freqs, Pxx, window_idx(t))
    name, freqs, Pxx = y_spectra_by_name[name]
    y_fundamental = get_max_frequency(freqs, Pxx, window_idx(t))

    name, freqs, Pxx = z_spectra_by_name[name]
    z_fundamental = get_max_frequency(freqs, Pxx, window_idx(t))

    try:
        x = butter_bandpass_filter(ts['x'], x_fundamental - space, x_fundamental + space, Fs)
        y = butter_bandpass_filter(ts['y'], y_fundamental - space, y_fundamental + space, Fs)
        z = butter_bandpass_filter(ts['z'], z_fundamental - space, z_fundamental + space, Fs)
    except ValueError:
        continue

    ax = plt.figure().add_subplot(projection='3d')
    ax.scatter(x, y, z, zdir='z', s=1, color='navy')
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    ax.set_zlabel('z')
    ax.set_xlim(-2, 2)
    ax.set_ylim(-2, 2)
    ax.set_zlim(-2, 2)
    ax.zaxis.labelpad = -0.7
    plt.show()

	x	y	z
count	250000	250000	250000
mean	0.00442377	0.00256424	0.00729225
std	0.628605	0.0284248	0.409702
min	-2.48307	-0.0753402	-1.50922
25%	-0.408632	-0.0203851	-0.279686
50%	0.0140406	0.00369867	-0.00664152
75%	0.434552	0.0257897	0.281814
max	2.06348	0.12156	2.09294

	x	y	z
count	250000	250000	250000
mean	0.0091607	0.00282899	0.00110133
std	0.108023	0.0304018	0.481472
min	-0.453353	-0.108466	-1.76433
25%	-0.0582322	-0.0188245	-0.336619
50%	0.0179518	0.00321832	-0.0209285
75%	0.0850128	0.025359	0.31991
max	0.322083	0.0957236	2.19127

	x	y	z
count	250000	250000	250000
mean	0.00812089	0.0028024	0.0175564
std	0.199627	0.0420892	1.9319
min	-0.761782	-0.147555	-4.2895
25%	-0.125994	-0.0262874	-1.70635
50%	0.0210197	0.00320514	-0.0998429
75%	0.155579	0.0328518	1.77266
max	0.513142	0.136843	4.5501

	x	y	z
count	250000	250000	250000
mean	0.0782755	0.00476235	0.0240576
std	1.57258	0.489683	1.94964
min	-3.40253	-1.30722	-3.30159
25%	-1.20947	-0.344713	-1.82606
50%	-0.0172133	0.0869372	-0.178371
75%	1.32615	0.470354	1.81411
max	4.26328	0.649884	4.78323

	x	y	z
count	250000	250000	250000
mean	0.0302217	0.0155812	-0.000449872
std	1.5222	2.99196	1.80557
min	-4.78097	-9.09082	-7.13457
25%	-1.122	-2.14442	-1.26887
50%	0.130627	0.558218	-0.052664
75%	1.0827	2.38637	1.39874
max	4.67537	5.6701	4.50713

	f	y	1x
0	253.296	0.178571	1
1	3.05176	0.0797164	0.0120482
2	274.658	0.0583556	1.08434
3	360.107	0.0343078	1.42169
4	555.42	0.0267719	2.19277

	f	y	1x
0	9.15527	0.0321657	1
1	125.122	0.0114684	13.6667
2	18.3105	0.0070992	2
3	708.008	0.00466999	77.3333
4	466.919	0.00390644	51

	f	y	1x
0	42.7246	0.104049	1
1	253.296	0.0174547	5.92857
2	338.745	0.0091005	7.92857
3	234.985	0.00612845	5.5
4	158.691	0.00422818	3.71429

	f	y	1x
0	42.7246	0.829558	1
1	3.05176	0.199422	0.0714286
2	115.967	0.0491828	2.71429
3	292.969	0.025542	6.85714
4	546.265	0.0137255	12.7857

	f	y	1x
0	347.9	0.647856	1
1	207.52	0.191531	0.596491
2	137.329	0.139324	0.394737
3	488.281	0.0564564	1.40351
4	357.056	0.0440482	1.02632

	f	y	1x
0	42.7246	0.447526	1
1	3.05176	0.116225	0.0714286
2	115.967	0.0666007	2.71429
3	198.364	0.0373502	4.64286
4	631.714	0.0242998	14.7857

	f	y	1x
0	42.7246	0.00383932	1
1	97.6562	0.00269748	2.28571
2	61.0352	0.00229201	1.42857
3	466.919	0.00211743	10.9286
4	73.2422	0.00181217	1.71429

	f	y	1x
0	9.15527	0.00710433	1
1	122.07	0.00281931	13.3333
2	94.6045	0.0021457	10.3333
3	143.433	0.00198247	15.6667
4	18.3105	0.00187833	2

	f	y	1x
0	42.7246	0.0171923	1
1	210.571	0.00400206	4.92857
2	234.985	0.00266793	5.5
3	299.072	0.00210059	7
4	338.745	0.0018251	7.92857

	f	y	1x
0	42.7246	0.34031	1
1	210.571	0.0167542	4.92857
2	421.143	0.0108792	9.85714
3	61.0352	0.00648204	1.42857
4	183.105	0.00478527	4.28571

	x	y	z
count	250000	250000	250000
mean	0.0438275	0.00369081	0.0156163
std	1.05157	0.253477	1.90162
min	-2.86968	-0.953974	-4.21464
25%	-0.725317	-0.164062	-1.71735
50%	-0.0285322	0.0863153	-0.00410187
75%	0.779466	0.21807	1.68526
max	3.5737	0.311572	4.68837

	f	y	1x
0	347.9	1.13694	1
1	70.1904	0.249652	0.201754
2	213.623	0.144129	0.614035
3	106.812	0.0789232	0.307018
4	256.348	0.0629811	0.736842

	f	y	1x
0	42.7246	0.14479	1
1	256.348	0.0180443	6
2	12.207	0.00889324	0.285714
3	344.849	0.00475763	8.07143
4	430.298	0.00405999	10.0714

	f	y	1x
0	42.7246	0.140939	1
1	2706.91	0.0292876	63.3571
2	424.194	0.0229368	9.92857
3	115.967	0.0209491	2.71429
4	2496.34	0.0165197	58.4286

	f	y	1x
0	42.7246	0.129175	1
1	125.122	0.0385745	2.92857
2	3421.02	0.0225597	80.0714
3	3469.85	0.0203046	81.2143
4	2966.31	0.016839	69.4286

	f	y	1x
0	42.7246	1.31154	1
1	6.10352	0.0558989	0.142857
2	515.747	0.0288995	12.0714
3	253.296	0.0202677	5.92857
4	2771	0.0182706	64.8571

	f	y	1x
0	42.7246	1.27444	1
1	503.54	0.0957012	11.7857
2	1007.08	0.0404181	23.5714
3	924.683	0.0261103	21.6429
4	1049.8	0.0202111	24.5714

	f	y	1x
0	36.6211	0.736213	1
1	1080.32	0.10298	29.5
2	247.192	0.0780015	6.75
3	1116.94	0.0704766	30.5
4	70.1904	0.0567678	1.91667

	f	y	1x
0	42.7246	1.25633	1
1	503.54	0.0663159	11.7857
2	967.407	0.0507056	22.6429
3	546.265	0.0313838	12.7857
4	1000.98	0.0270391	23.4286

Measurement exploratory data analysis¶